!##############################################################################
! MODULE globals
!
! This code defines GLOBALS AND FUNCTIONS to solve life-cycle problem 
! with couples, singles, and portfolio choice
!
! Author: Annika Bacher (annika.bacher@bi.no)   
!
! This version: May 2024
!
!##############################################################################

include "toolbox_own.f90"

module globals

        use toolbox_own

        implicit none

        
        !  ********  gender specific ********

        ! 1) model structural parameters 
        
        real*8 :: betaf, betam, betac, p_rec
        real*8 :: gammaf, gammam, gammac

        integer, parameter :: JJ = 36, minage = 29
        real*8, parameter  :: deltas = 0.73d0, deltac = 0.73d0 ! ROBUSTNESS: change to 73d0  
        real*8, parameter  :: tau = 0.27d0  ! (Kaplan & Huggett 2016)
        
        ! 2) SMP costs 

        real*8 :: SMFc(JJ), SMFs(JJ)
        
        real*8 :: step, SF_young, SF_old 
        
        integer :: numb
        
        ! 3) do-loops 
        
        integer :: j, zi, ti, l, i, k
        
        ! do-loops for hc return (variables prime)
        
        integer :: zpr, rpr, apr 
        
        ! 4) female income process by marital status (s: single, c: couple)
        
        ! determinstic part

        integer, parameter :: naa = 2
        real*8 :: dist_thetaf(naa), thetafs(naa), efffs(JJ+1),yfs, eff1fs, eff2fs, & 
                                thetafc(naa), efffc(JJ+1),yfc, eff1fc, eff2fc
        

        ! stochastic part, persistent component

        integer, parameter :: nz = 3
		
		real*8 :: rhofs, zfs_boom(nz), zfs_rec(nz), varfs_boom, varfs_rec, mufs_boom, mufs_rec
	
        real*8 :: rhofc, zfc_boom(nz), zfc_rec(nz), varfc_boom, varfc_rec, mufc_boom, mufc_rec

        ! stochastic part, transitory component

        real*8 :: efs(nz), rhoefs, varefs, muefs, rhoefc, varefc, muefc

        ! combining both shocks

        integer, parameter :: nzz = nz*nz

        real*8 :: afs_boom(nzz), afs_rec(nzz)
        
        real*8, dimension(nzz,2) :: afs_comb

        real*8, dimension(nz, nz) :: pizfs_boom, pizfs_rec, piefs, pizfc_boom, pizfc_rec

        real*8, dimension(nzz, nzz) ::  pi1fs_boom,pi1fs_rec, pi2fs, pifs_boom, pifs_rec                           
                                        
        
        ! 5) male income process

        ! deterministic part

        real*8 :: dist_thetam(naa), thetams(naa), effms(JJ+1), yms, eff1ms, eff2ms, & 
                             thetamc(naa), effmc(JJ+1), ymc, eff1mc, eff2mc
        
        ! stochastic part, persistent component

		real*8 :: rhoms, zms_boom(nz), zms_rec(nz), varms_boom, varms_rec, mums_boom, mums_rec

        real*8 :: rhomc, zmc_boom(nz), zmc_rec(nz), varmc_boom, varmc_rec, mumc_boom, mumc_rec

        ! stochastic part, transitory component

        real*8 :: ems(nz), rhoems, varems, muems, rhoemc, varemc, muemc

        ! combining both shocks

        real*8 :: ams_boom(nzz), ams_rec(nzz)
        
        real*8, dimension(nzz,2) :: ams_comb

        real*8, dimension(nz, nz)   ::  pizms_boom, pizms_rec, piems, pizmc_boom, pizmc_rec

        real*8, dimension(nzz, nzz) ::  pi1ms_boom, pi1ms_rec, pi2ms, pims_boom, pims_rec
        
        
        ! 6) joint income process for couples
        
        real*8, dimension(nz*nz,2) :: coupgrid
        
        real*8, dimension(nz*nz, nz*nz) :: Pcoup, picpers
        
        real*8, dimension(nzz*nzz, 2) :: coupgrid_all_boom, coupgrid_all_rec
        
        integer, dimension(2) :: nn = (/nz, nz/)
        
        real*8 :: am_c(nzz), af_c(nzz), corr, sigmac(2), mucoup(2)
        
        real*8, dimension(nzz*nzz,nzz*nzz) :: pic_boom, pic_rec, pi1c, pi2c_boom, pi2c_rec, pi3c_boom, pi3c_rec
        
        
        !  ********  not gender specific ********
                                     
        ! 1) asset returns
        
        real*8  :: r

        integer, parameter :: nrr = 4

        real*8 :: rhor, varr,mur_boom,mur_rec

        real*8, dimension(nrr,nrr) :: pi_risk_boom, pi_risk_rec
        
        real*8, dimension(nrr)     :: pir_boom, pir_rec
        
        real*8 :: ar_boom(nrr),ar_rec(nrr), ard_boom(nrr),ard_rec(nrr)
        
        real*8, dimension(nrr,2)   :: ard_comb
        
        ! 2) risky asset grid

        integer, parameter :: nrisk = 20   
        
        real*8 :: rgrid(nrisk), risk_min, risk_max
         
        ! 3) asset grid

        integer, parameter :: ncc = 100 

        real*8 :: cgrid(ncc), cmax, cmin, cstep, c_lb, c_ub
        
        real*8 :: add_low(JJ), add_high(JJ)
        
        ! 4) equivalence scales 
        
        real*8, dimension(JJ) :: nuf , num , nuc
        
        ! 5) divorce and marriage probability+ marriage market
        
        real*8, dimension(JJ,naa,nzz) :: mu_m, mu_f 
        
        real*8, dimension(JJ,naa,naa,nzz,nzz) :: lambda 
        
        real*8, dimension(naa*nzz,naa*nzz) :: trans_part, trans_part_m, trans_part_f
        
        real*8 :: mwealth(2,2,JJ+1)
                         
        ! 6) Value Functions and Optimization

        real*8, dimension(nrisk,ncc,nzz,naa,JJ)         :: Vopt_f , Vopt_m, cash_opt_m, cash_opt_f
                                                                                            
        real*8, dimension(ncc,nzz,naa,JJ)               :: VF_f, VF_m

        real*8, dimension(nrisk,ncc,nzz,nzz,naa,naa,JJ) :: Vopt_c, Vcoup_f, Vcoup_m ,cash_opt_c                                

        real*8, dimension(ncc,nzz,nzz,naa,naa,JJ)       :: VF_c, VF_cf, VF_cm
         
        ! 7) Policy Functions when single
        
        real*8, dimension(nrisk,ncc,nzz,naa,JJ) :: copt_f, copt_m, ctopt_sing_f, &
                                                   ctopt_sing_m, ctopt_coup_f, ctopt_coup_m
        
        real*8, dimension(ncc,nzz,naa,JJ) :: sharepol_f, sharepol_m, riskypol_f, riskypol_m, savepol_f,savepol_m, &
                                             cpol_f, cpol_m, &
                                             cprimepol_f,cprimepol_m, ct_sing_f, ct_sing_m, ct_coup_f, ct_coup_m
                                             
        ! 8) Policy Functions when couple  
                                             
        real*8, dimension(JJ,naa,naa,ncc,nzz,nzz) :: cpol_c, savepol_c, sharepol_c, riskypol_c, cprimepol_c    

        real*8, dimension(JJ,naa,naa,ncc,nzz,nzz,nrisk) ::  copt_c 
        
        
        ! 9) Displaying the results
        
        real*8 :: age(JJ)
        
        
        ! 10) Shock paths: productivity (male+female), asset return, marital status, aggregate state

        integer, parameter :: TT=JJ+1, SS=25000
        
        integer            :: s, r0f(SS), r0m(SS), z0f(SS), z0m(SS), z0c(2*SS), SS1, SS2, pers_shock
       
        real*8  :: ZC_zf(nzz), rr, aa, ZC_zm(nzz), ZC_mf(naa*nzz), ZC_mm(naa*nzz), ZC_zc(nzz*nzz), &
                   ZC_mc(nzz), ZC_help(nzz), SS1_h, SS2_h, ZC_rf(nrr), ZC_rm(nrr), ZC_st(nrr)
        
        integer, dimension(SS,TT) :: chain_retm, chain_retf, chain_prodf,chain_prodm, &
                                     chain_martf, chain_martm, chain_prodmf, chain_prodmm, chain_agg
    
        integer, dimension(2*SS,TT) :: chain_prodc
        
        ! 11) Simulation

        ! single men
        real*8, dimension(SS,TT) :: cash_path_sm, cprime_path_sm,&
                                    safe_path_sm, risky_path_sm, &
                                    share_path_sm, c_path_sm, SMP_path_sm, &
                                    cntsing_path_sm, cntcoup_path_sm, &
                                    hcval_path_sm, hcret_path_sm
                                                            
        ! single women                             
        real*8, dimension(SS,TT) :: cash_path_sf, cprime_path_sf, &
                                    safe_path_sf, risky_path_sf, &
                                    share_path_sf, c_path_sf, SMP_path_sf, &
                                    cntsing_path_sf, cntcoup_path_sf, &
                                    hcval_path_sf, hcret_path_sf
                                    
        ! couples (assets)
        real*8, dimension(2*SS,TT) :: cash_path_c, cprime_path_c, &
                                       safe_path_c, risky_path_c, share_path_c, c_path_c, SMP_path_c
                                   
        
        integer, dimension(SS,TT) :: gridl_path_m, gridh_path_m, grid_path_m, &
                                     gridl_path_f, gridh_path_f, grid_path_f
     
     
        integer, dimension(2*SS,TT) ::  grid_path_c
     
        real*8 :: AB2_f, AB2_m, draw, diff, inf, asset_sp, help
        
        integer :: AB_f, AB_m, asset_sp_ind 
        
        ! 12) cohort variables
        
        ! single men 
        real*8 :: cash_coh_sm(JJ), risky_coh_sm(JJ), safe_coh_sm(JJ), share_coh_sm(JJ), c_coh_sm(JJ), &
                                   SMP_coh_sm(JJ), ucshare_coh_sm(JJ), hcval_coh_sm(JJ), hcret_coh_sm(JJ)
                                   
        ! single women                      
        real*8 :: cash_coh_sf(JJ), risky_coh_sf(JJ), safe_coh_sf(JJ), share_coh_sf(JJ), c_coh_sf(JJ), &
                                   SMP_coh_sf(JJ), ucshare_coh_sf(JJ), hcval_coh_sf(JJ), hcret_coh_sf(JJ)
                                   
        ! couples (asset variables)                     
        real*8 :: cash_coh_c(JJ), risky_coh_c(JJ), safe_coh_c(JJ), &
                    share_coh_c(JJ), c_coh_c(JJ), SMP_coh_c(JJ), ucshare_coh_c(JJ)
            
                  
        ! 13) human capital returns
        
        real*8, dimension(JJ,nzz*ncc,nzz*ncc) :: dist_f_low, dist_f_high, dist_m_low, dist_m_high 
        
        real*8, dimension(JJ,ncc,nzz,JJ,ncc,nzz) :: sdf_f_low, sdf_f_high, sdf_m_low, sdf_m_high 
        
        real*8, dimension(JJ,ncc,nzz) :: hc_val_f_low, hc_val_f_high,hc_val_m_low, hc_val_m_high
        
                  
!##############################################################################  

!                           INCLUDED FUNCTIONS  
  
!##############################################################################  
  
        
contains

    ! 1) function for optimal intertemporal decision (single)

        function valfun_sing(cprime,share,gam,y,theta,eff,beta, &
                             c0,l,t,z,vfs,vc,pi_boom,pi_rec,mu,trans_part,nu,cont,SMF, ass_part, &
                             prod_boom,prod_rec,yi,yp,effi,effp,thetai,thetap,gend) result(valfres)
                             
        ! declare variables
            implicit none
    
            integer :: t, i, j, l, z, tpp, zpp, gend
    
            integer, dimension(nrr,nzz,2) :: clowst, chighst
            
            integer, dimension(nrr,2) :: clowst_ret, chighst_ret
            
            integer, dimension(nrr,nzz,nzz,naa,2)  ::  cash_mlow, cash_mhigh
            
            real*8  :: gg_s, gg_vc, valfres(4), gg_ret, ass_part(2)
    
            real*8  :: share,c0, cprime, consum, nu, SMF, coupstate_boom(2), coupstate_rec(2)
            
            real*8  :: y, yi, yp, theta, thetai, thetap(naa), eff, effi, effp, gam, beta, inf
    
            real*8  :: valf, mu, prod_boom(nzz), prod_rec(nzz), income_boom, income_rec
    
            real*8  :: exp_s, exp_ic, exp_ret, cont(ncc)
			
			real*8, dimension(nrr,nzz,2) :: cashtmrst
			
			real*8, dimension(nrr,2) :: cashtmrst_ret
    
            real*8, dimension(ncc,nzz) :: vfs
    
            real*8, dimension(nzz,nzz) :: pi_boom, pi_rec
    
            real*8, dimension(ncc,nzz,nzz,naa) :: vc
            
            real*8, dimension(naa*nzz,naa*nzz) :: trans_part
            
            real*8, dimension(nrr,nzz,nzz,naa,2) :: cash_m
            

        ! define infinity
            inf = 900000000000000000000000000d0

        ! HH chooses cprime
        ! cashtmr is assets cash-on-hand next period: assets + labor income
        
        if(l<JJ) then

            do i=1,nrr
				do j=1,nzz
    
            cashtmrst(i,j,1) = y*theta*eff*prod_boom(j)+ ((1d0-share)*(1d0+r)+share*ard_boom(i))*cprime
			cashtmrst(i,j,2) = y*theta*eff*prod_rec(j) + ((1d0-share)*(1d0+r)+share*ard_rec(i))*cprime
			
            clowst(i,j,1)    = max(cntarray(cgrid,ncc,cashtmrst(i,j,1)),1)
            chighst(i,j,1)   = min(clowst(i,j,1) + 1,ncc)
            if (clowst(i,j,1) == ncc) clowst(i,j,1)  = clowst(i,j,1)  -1 
            
            clowst(i,j,2)    = max(cntarray(cgrid,ncc,cashtmrst(i,j,2)),1)
            chighst(i,j,2)   = min(clowst(i,j,2) + 1,ncc)
            if (clowst(i,j,2) == ncc) clowst(i,j,2)  = clowst(i,j,2)  -1 

				end do
            end do
            
        ! linear interpolation of value function    
        ! specify the continuation values 
        
		! when staying single
            gg_s  = 0d0
            exp_s = 0d0
               
            do i = 1,nrr
				do j = 1, nzz
				
            gg_s = gg_s + ((1d0-p_rec)*pir_boom(i)*pi_boom(z,j)*(vfs(clowst(i,j,1),j) + &
                (cashtmrst(i,j,1)-cgrid(clowst(i,j,1)))*(vfs(chighst(i,j,1),j) - &
                vfs(clowst(i,j,1),j))/(cgrid(chighst(i,j,1)) - cgrid(clowst(i,j,1)))) + & 
                p_rec*pir_rec(i)*pi_rec(z,j)*(vfs(clowst(i,j,2),j) + &
                (cashtmrst(i,j,2)-cgrid(clowst(i,j,2)))*(vfs(chighst(i,j,2),j) - &
                vfs(clowst(i,j,2),j))/(cgrid(chighst(i,j,2)) - cgrid(clowst(i,j,2)))))    
                
                end do
            end do
            
            exp_s = gg_s
            valfres(3) = exp_s
            
            
            ! individual being in a couple 
            ! asset allocation when getting married: average assets, dep. on education
            
			do tpp = 1,naa		   ! partner education
				do j = 1, nzz      ! own productivity
				   do zpp = 1, nzz  ! partner productivity
			
			! determine income within couple
			if (gend ==1) then ! men 
			
			  if(zpp <= 3) then
              coupstate_boom(:) = coupgrid_all_boom(nz*(j-1)+zpp,:)
              coupstate_rec(:)  =  coupgrid_all_rec(nz*(j-1)+zpp,:)
              else if (zpp > 3 .AND. zpp < 7) then
              coupstate_boom(:) = coupgrid_all_boom((nz*(j-1)+zpp-nz)+nz*nzz,:)
              coupstate_rec(:)  =  coupgrid_all_rec((nz*(j-1)+zpp-nz)+nz*nzz,:)
              else if  (zpp > 6) then
              coupstate_boom(:) = coupgrid_all_boom((nz*(j-1)+zpp-2*nz)+2*nz*nzz,:)
              coupstate_rec(:)  =  coupgrid_all_rec((nz*(j-1)+zpp-2*nz)+2*nz*nzz,:)
              end if 
              			
		    income_boom   = yi*effi*thetai*coupstate_boom(2)+yp*effp*thetap(tpp)*coupstate_boom(1) 
		    income_rec    = yi*effi*thetai*coupstate_rec(2) +yp*effp*thetap(tpp)*coupstate_rec(1) 
		    
			elseif (gend ==2) then ! women
			
			  if(j <= 3) then
              coupstate_boom(:) = coupgrid_all_boom(nz*(zpp-1)+j,:)
              coupstate_rec(:)  =  coupgrid_all_rec(nz*(zpp-1)+j,:)
              else if (j > 3 .AND. j < 7) then
              coupstate_boom(:) = coupgrid_all_boom((nz*(zpp-1)+j-nz)+nz*nzz,:)
              coupstate_rec(:)  =  coupgrid_all_rec((nz*(zpp-1)+j-nz)+nz*nzz,:)
              else if  (j > 6) then
              coupstate_boom(:) = coupgrid_all_boom((nz*(zpp-1)+j-2*nz)+2*nz*nzz,:)
              coupstate_rec(:)  =  coupgrid_all_rec((nz*(zpp-1)+j-2*nz)+2*nz*nzz,:)
              end if 

			income_boom   = yi*effi*thetai*coupstate_boom(1)+yp*effp*thetap(tpp)*coupstate_boom(2) 
			income_rec    = yi*effi*thetai*coupstate_rec(1) +yp*effp*thetap(tpp)*coupstate_rec(2) 
			end if
			
			
			do i = 1, nrr
		
			cash_m(i,j,zpp,tpp,1)     = income_boom + ((1d0-share)*(1d0+r)+share*ard_boom(i))*cprime + ass_part(tpp) 		
			cash_m(i,j,zpp,tpp,2)     = income_rec  + ((1d0-share)*(1d0+r)+share*ard_rec(i))*cprime  + ass_part(tpp)
								   
            cash_mlow(i,j,zpp,tpp,1)   = max(cntarray(cgrid,ncc,cash_m(i,j,zpp,tpp,1)),1)
            cash_mhigh(i,j,zpp,tpp,1)  = min(cash_mlow(i,j,zpp,tpp,1) + 1,ncc)
            
            if (cash_mlow(i,j,zpp,tpp,1) == ncc) cash_mlow(i,j,zpp,tpp,1) = cash_mlow(i,j,zpp,tpp,1) -1
            
            cash_mlow(i,j,zpp,tpp,2)   = max(cntarray(cgrid,ncc,cash_m(i,j,zpp,tpp,2)),1)
            cash_mhigh(i,j,zpp,tpp,2)  = min(cash_mlow(i,j,zpp,tpp,2) + 1,ncc)
            
            if (cash_mlow(i,j,zpp,tpp,2) == ncc) cash_mlow(i,j,zpp,tpp,2) = cash_mlow(i,j,zpp,tpp,2) -1
            
            
			end do
			
						end do
                     end do 
               end do 
      
           
           ! specify continuation value, linear interpolation of value function
           gg_vc = 0d0  
           exp_ic = 0d0

                do tpp=1,naa
                   do i = 1,nrr
					  do j=1,nzz   			 ! own productivity (tomorrow)
						do zpp=1,nzz  	     ! partner productivity	

            gg_vc = gg_vc + ((1d0-p_rec)*pir_boom(i)*pi_boom(z,j)*trans_part((t-1)*nzz+z,(tpp-1)*nzz+zpp)* ( &
            vc(cash_mlow(i,j,zpp,tpp,1),j,zpp,tpp) + (cash_m(i,j,zpp,tpp,1)-cgrid(cash_mlow(i,j,zpp,tpp,1)))* &
            (vc(cash_mhigh(i,j,zpp,tpp,1),j,zpp,tpp) - vc(cash_mlow(i,j,zpp,tpp,1),j,zpp,tpp))/ &
            (cgrid(cash_mhigh(i,j,zpp,tpp,1)) - cgrid(cash_mlow(i,j,zpp,tpp,1)))) + &
						p_rec*pir_rec(i)*pi_rec(z,j)*trans_part((t-1)*nzz+z,(tpp-1)*nzz+zpp)* ( &
            vc(cash_mlow(i,j,zpp,tpp,2),j,zpp,tpp) + (cash_m(i,j,zpp,tpp,2)-cgrid(cash_mlow(i,j,zpp,tpp,2)))* &
            (vc(cash_mhigh(i,j,zpp,tpp,2),j,zpp,tpp) - vc(cash_mlow(i,j,zpp,tpp,2),j,zpp,tpp))/ &
            (cgrid(cash_mhigh(i,j,zpp,tpp,2)) - cgrid(cash_mlow(i,j,zpp,tpp,2)))) )

            
							end do
                        end do
                    end do
                end do
            
            exp_ic = gg_vc
            valfres(4) = exp_ic
   
                
            else    ! last period before retirment
            
            do i=1,nrr
    
            cashtmrst_ret(i,1) = y*theta*eff*prod_boom(z)*0.6d0 + ((1d0-share)*(1d0+r)+share*ard_boom(i))*cprime
			cashtmrst_ret(i,2) = y*theta*eff*prod_boom(z)*0.6d0 + ((1d0-share)*(1d0+r)+share*ard_rec(i))*cprime
			
            clowst_ret(i,1)    = max(cntarray(cgrid,ncc,cashtmrst_ret(i,1)),1)
            chighst_ret(i,1)   = min(clowst_ret(i,1) + 1,ncc)
            
            if (clowst_ret(i,1) == ncc) clowst_ret(i,1) = clowst_ret(i,1) -1
            
            clowst_ret(i,2)    = max(cntarray(cgrid,ncc,cashtmrst_ret(i,2)),1)
            chighst_ret(i,2)   = min(clowst_ret(i,2) + 1,ncc)

			if (clowst_ret(i,2) == ncc) clowst_ret(i,2) = clowst_ret(i,2) -1
            
            end do
            
            gg_ret  = 0d0
            exp_ret = 0d0
                
            do i = 1,nrr
            gg_ret = gg_ret+ ((1d0-p_rec)*pir_boom(i)*(cont(clowst_ret(i,1)) + &
                (cashtmrst_ret(i,1)-cgrid(clowst_ret(i,1)))*(cont(chighst_ret(i,1)) - &
                cont(clowst_ret(i,1)))/(cgrid(chighst_ret(i,1)) - cgrid(clowst_ret(i,1)))) + &
                p_rec*pir_rec(i)*(cont(clowst_ret(i,2)) + &
                (cashtmrst_ret(i,2)-cgrid(clowst_ret(i,2)))*(cont(chighst_ret(i,2)) - &
                cont(clowst_ret(i,2)))/(cgrid(chighst_ret(i,2)) - cgrid(clowst_ret(i,2)))))
            end do
            
            exp_ret=gg_ret
                        
            valfres(3) = exp_ret
            valfres(4) = exp_ret
            
            end if
            
            
        ! compute consumption
            if(share==0d0) then
            consum = c0-cprime
            else
            consum = c0-SMF-cprime
            end if
            
        
		! robustness: receving bequest at age 55 (5000 is normalization value)

	    !if(l==26 .AND. gend == 1) consum = consum + 1.83d0
	    !if(l==26 .AND. gend == 2) consum = consum + 1.60d0

            
        ! impossible allocations ('consumption floor')
            if(consum<=0.0002d0) consum = 0.0002d0

         
        ! compute the value function 
        
            if(l<JJ) then
  
            ! standard CRRA 
                                        
            valf = nu*((consum/nu)**(1d0-gam))/(1d0-gam) +beta*(1d0-mu)*exp_s + &
                                    beta*mu*exp_ic
            
            else
            
            ! period before retirement

            valf = nu*((consum/nu)**(1d0-gam))/(1d0-gam) + beta*exp_ret
            end if
                        
                        
            valfres(1) = valf
            valfres(2) = consum
            
            end function valfun_sing
            
    ! 4) function for optimal intertemporal decision (couple)

        function valfun_coup(cprime,share,gami,gamp,yi,yp,thetai,thetap,effi,effp, &
                             betai,betap,c0,l,zi,zp,vfi,vfp,vfc,vci,vcp, &
                             pii_boom,pii_rec,pip_boom,pip_rec,pic_boom,pic_rec,lambda,nu, &
                             cont,conti, contp,SMF,yf,ym,efff,effm,thetaf,thetam, &
                             prodf_boom, prodf_rec,prodm_boom,prodm_rec) result(val_results)
                                
        
        ! declare variables
        
            implicit none
            
            integer ::  zi,zp, i, l, k, j, m, n
            
            integer, dimension(nrr,2) :: clowst_ret, chighst_ret
    
            integer, dimension(nrr,nzz,nzz,2) ::  clowst, chighst 
            
            integer, dimension(nrr,nzz,2) :: clowst_m, chighst_m, clowst_f, chighst_f

            real*8, dimension(nrr,nzz,nzz,2) :: cashtmrst
            
            real*8, dimension(nrr,nzz,2) :: cashtmrst_m,  cashtmrst_f 
            
            real*8, dimension(nrr,2) :: cashtmrst_ret
            
            real*8  ::  gg_ret, gg_reti, gg_retp 
            
            real*8  ::  gg_si, gg_sp, val_results(4), coupstate_boom(2), coupstate_rec(2)
    
            real*8  ::  share, c0, cprime, cprime_m, cprime_f, consum, val_ci, val_cp, inf, SMF
            
            real*8  ::  yi, yp, thetai, thetap, effi, effp, gami, gamp, betai, betap, nu
            
            real*8  ::  yf, ym, efff, effm, thetaf, thetam 
            
            real*8  ::  prodf_boom(nzz), prodf_rec(nzz), prodm_boom(nzz), prodm_rec(nzz)
    
            real*8  ::  valf, gg_c, gg_vcp, gg_vci
            
            real*8  ::  exp_si, exp_sp ,exp_c, exp_vci, exp_vcp, lambda, exp_ret, exp_reti, exp_retp
    
            real*8, dimension(ncc,nzz) :: vfi, vfp
            
            real*8, dimension(ncc,nzz,nzz) ::  vfc, vci, vcp
            
            real*8, dimension(ncc) :: cont, conti, contp
    
            real*8, dimension(nzz,nzz) ::  pii_boom, pii_rec, pip_boom, pip_rec
            
            real*8, dimension(nzz*nzz,nzz*nzz) :: pic_boom, pic_rec
            
         
        ! define infinity
            inf = 900000000000000000000000000d0
        
        ! HH chooses cprime, specify cash tomporrow (dep. on states)
        
        if(l<JJ) then
        
        ! in case of seperation, assets are split 50-50 (with 10% destroyed)
            cprime_m = cprime*0.45d0  !cprime*0.315d0 ! ROBUSTNESS, MODIFIED ASSET SPLITTING RULE
            cprime_f = cprime*0.45d0  !cprime*0.585d0

    
		! case of staying married
            do i=1,nrr
				do j=1,nzz     ! productivity of i (woman)
					do m=1,nzz ! productivity of p (man)
					
		 !specify individual productivity states of spouses (j refers to wife, m to husband)     
              if(j <= 3) then
              coupstate_boom(:) = coupgrid_all_boom(nz*(m-1)+j,:)
              coupstate_rec(:)  =  coupgrid_all_rec(nz*(m-1)+j,:)
              else if (j > 3 .AND. j < 7) then
              coupstate_boom(:) = coupgrid_all_boom((nz*(m-1)+j-nz)+nz*nzz,:)
              coupstate_rec(:)  =  coupgrid_all_rec((nz*(m-1)+j-nz)+nz*nzz,:)
              else if  (j > 6) then
              coupstate_boom(:) = coupgrid_all_boom((nz*(m-1)+j-2*nz)+2*nz*nzz,:)
              coupstate_rec(:)  =  coupgrid_all_rec((nz*(m-1)+j-2*nz)+2*nz*nzz,:)
              end if 
    
            cashtmrst(i,j,m,1) = yi*thetai*effi*coupstate_boom(1)+yp*thetap*effp*coupstate_boom(2) &
							   + ((1d0-share)*(1d0+r)+share*ard_boom(i))*cprime
							   
		    cashtmrst(i,j,m,2) = yi*thetai*effi*coupstate_rec(1)+yp*thetap*effp*coupstate_rec(2) &
							   + ((1d0-share)*(1d0+r)+share*ard_rec(i))*cprime
    
            clowst(i,j,m,1)    = max(cntarray(cgrid,ncc,cashtmrst(i,j,m,1)),1)
            chighst(i,j,m,1)   = min(clowst(i,j,m,1) + 1,ncc)
            
            if (clowst(i,j,m,1) == ncc) clowst(i,j,m,1)  = clowst(i,j,m,1)  -1 
            
            clowst(i,j,m,2)    = max(cntarray(cgrid,ncc,cashtmrst(i,j,m,2)),1)
            chighst(i,j,m,2)   = min(clowst(i,j,m,2) + 1,ncc)
            
            if (clowst(i,j,m,2) == ncc) clowst(i,j,m,2)  = clowst(i,j,m,2)  -1
            
					end do
				end do
            end do
            
            
            ! case of separation
            do i=1,nrr
				do j=1,nzz
    
            cashtmrst_m(i,j,1) = ym*effm*thetam*prodm_boom(j)+ ((1d0-share)*(1d0+r)+share*ard_boom(i))*cprime_m
            cashtmrst_m(i,j,2) = ym*effm*thetam*prodm_rec(j) + ((1d0-share)*(1d0+r)+share*ard_rec(i))*cprime_m
            
            cashtmrst_f(i,j,1) = yf*efff*thetaf*prodf_boom(j)+ ((1d0-share)*(1d0+r)+share*ard_boom(i))*cprime_f
            cashtmrst_f(i,j,2) = yf*efff*thetaf*prodf_rec(j) + ((1d0-share)*(1d0+r)+share*ard_rec(i))*cprime_f
    
            clowst_m(i,j,1)    = max(cntarray(cgrid,ncc,cashtmrst_m(i,j,1)),1)
            clowst_m(i,j,2)    = max(cntarray(cgrid,ncc,cashtmrst_m(i,j,2)),1)
            
            clowst_f(i,j,1)    = max(cntarray(cgrid,ncc,cashtmrst_f(i,j,1)),1)
            clowst_f(i,j,2)    = max(cntarray(cgrid,ncc,cashtmrst_f(i,j,2)),1)
            
            chighst_m(i,j,1)   = min(clowst_m(i,j,1) + 1,ncc)
            chighst_m(i,j,2)   = min(clowst_m(i,j,2) + 1,ncc)
            
            chighst_f(i,j,1)   = min(clowst_f(i,j,1) + 1,ncc)
            chighst_f(i,j,2)   = min(clowst_f(i,j,2) + 1,ncc)
            
            if (clowst_m(i,j,1) == ncc) clowst_m(i,j,1)  = clowst_m(i,j,1)  -1 
            if (clowst_m(i,j,2) == ncc) clowst_m(i,j,2)  = clowst_m(i,j,2)  -1 
           
            if (clowst_f(i,j,1) == ncc) clowst_f(i,j,1)  = clowst_f(i,j,1)  -1
			if (clowst_f(i,j,2) == ncc) clowst_f(i,j,2)  = clowst_f(i,j,2)  -1
			
				end do
            end do


        ! continuation values: linear interpolation of value function

        ! Value of single (seperately for individual i and partner p)
			gg_si  = 0d0
			exp_si = 0d0
			
			gg_sp  = 0d0
			exp_sp = 0d0

            do i = 1,nrr
				do j = 1,nzz
            gg_si = gg_si+ ( (1d0-p_rec)*pir_boom(i)*pii_boom(zi,j)*(vfi(clowst_f(i,j,1),j) + &
               (cashtmrst_f(i,j,1)-cgrid(clowst_f(i,j,1)))*(vfi(chighst_f(i,j,1),j) - &
               vfi(clowst_f(i,j,1),j))/(cgrid(chighst_f(i,j,1)) - cgrid(clowst_f(i,j,1)))) + &
               p_rec*pir_rec(i)*pii_rec(zi,j)*(vfi(clowst_f(i,j,2),j) + &
               (cashtmrst_f(i,j,2)-cgrid(clowst_f(i,j,2)))*(vfi(chighst_f(i,j,2),j) - &
               vfi(clowst_f(i,j,2),j))/(cgrid(chighst_f(i,j,2)) - cgrid(clowst_f(i,j,2)))))
               end do
            end do 
         
            do i = 1,nrr
				do j = 1,nzz
            gg_sp = gg_sp+ ( (1d0-p_rec)*pir_boom(i)*pip_boom(zp,j)*(vfp(clowst_m(i,j,1),j) + &
               (cashtmrst_m(i,j,1)-cgrid(clowst_m(i,j,1)))*(vfp(chighst_m(i,j,1),j) - &
               vfp(clowst_m(i,j,1),j))/(cgrid(chighst_m(i,j,1)) - cgrid(clowst_m(i,j,1)))) + &
               p_rec*pir_rec(i)*pip_rec(zp,j)*(vfp(clowst_m(i,j,2),j) + &
               (cashtmrst_m(i,j,2)-cgrid(clowst_m(i,j,2)))*(vfp(chighst_m(i,j,2),j) - &
               vfp(clowst_m(i,j,2),j))/(cgrid(chighst_m(i,j,2)) - cgrid(clowst_m(i,j,2)))) )
               end do
            end do 
        
            exp_si=gg_si
			exp_sp=gg_sp
  
		   ! Value for individual of being in a couple (seperately for i and p)
			gg_vci  = 0d0
			exp_vci = 0d0
			
			gg_vcp  = 0d0
			exp_vcp = 0d0
			   
           ! first: position on joint productivity grid TODAY, given individual states
           if(zi <= 3) then
           k = nz*(zp-1)+zi
           else if (zi > 3 .AND. zi < 7) then
           k = (nz*(zp-1)+zi-nz)+nz*nzz
           else if  (zi > 6) then
           k = (nz*(zp-1)+zi-2*nz)+2*nz*nzz
           end if 
           
           ! second: "translate" back into two distinct state variables
           do m=1,nz 
             do j=1,nzz
               do i= (m-1)*nzz*nz+(j-1)*nz+1,(m-1)*nzz*nz+nz*j
               
              do n = 1,nrr
              gg_vci = gg_vci+ ( (1d0-p_rec)*pir_boom(n)*pic_boom(k,i)* &
               (vci(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1),i-((m-1)*nzz*nz)-nz*(j-m),j) + &
               (cashtmrst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)-cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)))* &
               (vci(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1),i-((m-1)*nzz*nz)-nz*(j-m),j) - &
                vci(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1),i-((m-1)*nzz*nz)-nz*(j-m),j))/ &
               (cgrid(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)) - cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)))) + &
               p_rec*pir_rec(n)*pic_rec(k,i)* &
               (vci(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2),i-((m-1)*nzz*nz)-nz*(j-m),j) + &
               (cashtmrst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)-cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)))* &
               (vci(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2),i-((m-1)*nzz*nz)-nz*(j-m),j) - &
                vci(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2),i-((m-1)*nzz*nz)-nz*(j-m),j))/ &
               (cgrid(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)) - cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)))) )
              end do 
            

				end do
              end do   
            end do
            
            exp_vci = gg_vci
            
                        
            do m=1,nz
              do j=1,nzz
                do i= (m-1)*nzz*nz+(j-1)*nz+1,(m-1)*nzz*nz+nz*j
            
                 do n = 1,nrr
              gg_vcp = gg_vcp + ( (1d0-p_rec)*pir_boom(n)*pic_boom(k,i)* &
               (vcp(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),1),j,i-((m-1)*nzz*nz)-nz*(j-m)) + &
               (cashtmrst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),1)-cgrid(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),1)))* &
               (vcp(chighst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),1),j,i-((m-1)*nzz*nz)-nz*(j-m)) - &
                vcp(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),1),j,i-((m-1)*nzz*nz)-nz*(j-m)))/ &
               (cgrid(chighst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),1)) - cgrid(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),1)))) + &
               p_rec*pir_rec(n)*pic_rec(k,i)* &
               (vcp(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),2),j,i-((m-1)*nzz*nz)-nz*(j-m)) + &
               (cashtmrst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),2)-cgrid(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),2)))* &
               (vcp(chighst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),2),j,i-((m-1)*nzz*nz)-nz*(j-m)) - &
                vcp(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),2),j,i-((m-1)*nzz*nz)-nz*(j-m)))/ &
               (cgrid(chighst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),2)) - cgrid(clowst(n,j,i-((m-1)*nzz*nz)-nz*(j-m),2)))) )
                 end do 

                end do
              end do  
            end do
            
            exp_vcp = gg_vcp
            
		  ! value of couple
           gg_c   = 0d0
		   exp_c  = 0d0
		   
           ! first: position on joint productivity grid TODAY, given individual states
           if(zi <= 3) then
           k = nz*(zp-1)+zi
           else if (zi > 3 .AND. zi < 7) then
           k = (nz*(zp-1)+zi-nz)+nz*nzz
           else if  (zi > 6) then
           k = (nz*(zp-1)+zi-2*nz)+2*nz*nzz
           end if 
           
           ! second: "translate" back into two distinct state variables
           do m=1,nz
             do j=1,nzz
               do i= (m-1)*nzz*nz+(j-1)*nz+1,(m-1)*nzz*nz+nz*j
        
            do n = 1,nrr 
            gg_c= gg_c+ ( (1d0-p_rec)*pir_boom(n)*pic_boom(k,i)* &
			    (vfc(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1),i-((m-1)*nzz*nz)-nz*(j-m),j) + &
                (cashtmrst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)- &
                 cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)))* &
                (vfc(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1),i-((m-1)*nzz*nz)-nz*(j-m),j) - &
                 vfc(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1),i-((m-1)*nzz*nz)-nz*(j-m),j))/ &
                (cgrid(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)) - cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,1)))) + &
                p_rec*pir_rec(n)*pic_rec(k,i)* &
			    (vfc(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2),i-((m-1)*nzz*nz)-nz*(j-m),j) + &
                (cashtmrst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)- &
                 cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)))* &
                (vfc(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2),i-((m-1)*nzz*nz)-nz*(j-m),j) - &
                 vfc(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2),i-((m-1)*nzz*nz)-nz*(j-m),j))/ &
                (cgrid(chighst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)) - cgrid(clowst(n,i-((m-1)*nzz*nz)-nz*(j-m),j,2)))) )
            end do
            
				   end do
                end do
            end do
            
            exp_c = gg_c
         

            else       ! cont value during retirement 
            
            do i=1,nrr
    
            cashtmrst_ret(i,1) = yp*thetap*effp*coupgrid_all_boom(((zp-1)*nz) + (nzz*nz)*1+ 2,2)*0.6d0*1.7d0 &
							   + ((1d0-share)*(1d0+r)+share*ard_boom(i))*cprime
							   
		    cashtmrst_ret(i,2) = yp*thetap*effp*coupgrid_all_boom(((zp-1)*nz) + (nzz*nz)*1+ 2,2)*0.6d0*1.7d0 &
							   + ((1d0-share)*(1d0+r)+share*ard_rec(i))*cprime
    
            clowst_ret(i,1)    = max(cntarray(cgrid,ncc,cashtmrst_ret(i,1)),1)
            chighst_ret(i,1)   = min(clowst_ret(i,1) + 1,ncc)
            
            clowst_ret(i,2)    = max(cntarray(cgrid,ncc,cashtmrst_ret(i,2)),1)
            chighst_ret(i,2)   = min(clowst_ret(i,2) + 1,ncc)
            
            if (clowst_ret(i,1) == ncc) clowst_ret(i,1)  = clowst_ret(i,1)  -1 
            if (clowst_ret(i,2) == ncc) clowst_ret(i,2)  = clowst_ret(i,2)  -1
            
			end do

			! Value of couple during retirment
            gg_ret  = 0d0
            exp_ret = 0d0

            do i = 1,nrr 
            gg_ret = gg_ret + ((1d0-p_rec)*pir_boom(i)*(cont(clowst_ret(i,1)) + &
                (cashtmrst_ret(i,1)-cgrid(clowst_ret(i,1)))*(cont(chighst_ret(i,1)) - &
                cont(clowst_ret(i,1)))/(cgrid(chighst_ret(i,1)) - cgrid(clowst_ret(i,1))))+&
                p_rec*pir_rec(i)*(cont(clowst_ret(i,2)) + &
                (cashtmrst_ret(i,2)-cgrid(clowst_ret(i,2)))*(cont(chighst_ret(i,2)) - &
                cont(clowst_ret(i,2)))/(cgrid(chighst_ret(i,2)) - cgrid(clowst_ret(i,2)))) )
            end do
            
            exp_ret=gg_ret
            
            ! Value for individual of being in a couple during ret (i and p)
            gg_reti = 0d0
            gg_retp = 0d0
            
            exp_reti = 0d0
            exp_retp = 0d0
        
            do i = 1,nrr 
            gg_reti = gg_reti+ ((1d0-p_rec)*pir_boom(i)*(conti(clowst_ret(i,1)) + &
                (cashtmrst_ret(i,1)-cgrid(clowst_ret(i,1)))*(conti(chighst_ret(i,1)) - &
                conti(clowst_ret(i,1)))/(cgrid(chighst_ret(i,1)) - cgrid(clowst_ret(i,1)))) + &
                p_rec*pir_rec(i)*(conti(clowst_ret(i,2)) + &
                (cashtmrst_ret(i,2)-cgrid(clowst_ret(i,2)))*(conti(chighst_ret(i,2)) - &
                conti(clowst_ret(i,2)))/(cgrid(chighst_ret(i,2)) - cgrid(clowst_ret(i,2)))) )
            end do
        
            do i = 1,nrr 
            gg_retp = gg_retp+ ((1d0-p_rec)*pir_boom(i)*(contp(clowst_ret(i,1)) + &
                (cashtmrst_ret(i,1)-cgrid(clowst_ret(i,1)))*(contp(chighst_ret(i,1)) - &
                contp(clowst_ret(i,1)))/(cgrid(chighst_ret(i,1)) - cgrid(clowst_ret(i,1)))) +&
                p_rec*pir_rec(i)*(contp(clowst_ret(i,2)) + &
                (cashtmrst_ret(i,2)-cgrid(clowst_ret(i,2)))*(contp(chighst_ret(i,2)) - &
                contp(clowst_ret(i,2)))/(cgrid(chighst_ret(i,2)) - cgrid(clowst_ret(i,2)))) )
            end do
            
            exp_reti =gg_reti
            exp_retp =gg_retp
            end if

        !  compute consumption
            if(share==0d0) then
            consum = c0-cprime
            else
            consum = c0-SMF-cprime
            end if
        
        ! robustness: receving bequest at age 55 (50000 is normalization value)
	    !if(l==26) consum = consum + 3.40d0
     
        ! impossible allocations (conumption floor')
            if(consum<=0.0002d0) consum = 0.0002d0
 
        ! compute the value function and value of each individual of being in a couple
    
            if(l==JJ) then

            valf    =  nu*((consum/nu)**(1d0-gami))/(1d0-gami) + &
                            betai*exp_ret
            
            val_ci  =  nu*((consum/nu)**(1d0-gami))/(1d0-gami) + &
                            betai*exp_reti
                            
            val_cp  =  nu*((consum/nu)**(1d0-gamp))/(1d0-gamp) + &
                            betap*exp_retp
            
            
            else 
            
            
            valf    =  nu*((consum/nu)**(1d0-gamp))/(1d0-gamp) + &
                                                   betai*(1d0-lambda)*exp_c    + betai*lambda*(exp_si +exp_sp)
            
            val_ci  = nu*((consum/nu)**(1d0-gami))/(1d0-gami) + &
                                                    betai*(1d0-lambda)*exp_vci + betai*lambda*exp_si
                                                    
            val_cp  = nu*((consum/nu)**(1d0-gamp))/(1d0-gamp) + & 
                                                    betap*(1d0-lambda)*exp_vcp + betap*lambda*exp_sp
                                                    
            
                                                                                            
            end if  
                                
            val_results(1) = valf
            val_results(2) = val_ci
            val_results(3) = val_cp
            val_results(4) = consum

            
        end function valfun_coup
        
    ! 5) golden function (single) 
    
        function golden_sing(aa,bb,share,gam,y,theta,eff,beta, &
                        c0,l,t,z,vfs,vc,pi_boom,pi_rec,marr,trans_part,nu,cont,SMF, ass_part,prod_boom, prod_rec, &
                        yi,yp,effi,effp,thetai,thetap,gend) result(output)
                                

            ! declare variables
            implicit none
    
            real*8  :: alpha1, alpha2, d, x1, x2, tol, output(5), f1, f2
            
            real*8  :: aa, bb, prod_boom(nzz), prod_rec(nzz)
            
            integer :: t, l, z, gend
    
            real*8  :: share,c0,marr, valfres(4), ass_part(2)
            
            real*8  :: y, theta, eff, gam,  beta, nu, SMF 
            
            real*8  :: yi, yp, effi, effp, thetai, thetap(naa)

            real*8, dimension(ncc,nzz) :: vfs
            
            real*8, dimension(ncc) :: cont
            
            real*8, dimension(ncc,nzz,nzz,naa) :: vc
    
            real*8, dimension(nzz,nzz) :: pi_boom, pi_rec
            
            real*8, dimension(nzz*naa,nzz*naa) :: trans_part
            
            ! find optimal value with golden section search
    
            tol     =   0.000000001d0
            alpha1  =   (3d0-sqrt(5d0))/2d0
            alpha2  =   (sqrt(5d0)-1d0)/2d0
                d   =   bb-aa
                x1  =   aa+alpha1*d   
                x2  =   aa+alpha2*d   
                valfres  =   &
                valfun_sing(x1,share,gam,y,theta,eff,beta, &
                            c0,l,t,z,vfs,vc,pi_boom,pi_rec,marr,trans_part,nu,cont, SMF, ass_part, &
                            prod_boom, prod_rec,yi,yp,effi,effp,thetai,thetap,gend)
                                   
                f1  = valfres(1)
                
                if(f1==-900000000000000000000000000d0 .AND. share == 0d0) write(*,*) 'vf is infinity (for singles) f1'
                
                valfres  =   &
                valfun_sing(x2,share,gam,y,theta,eff,beta,&
                            c0,l,t,z,vfs,vc,pi_boom,pi_rec,marr,trans_part,nu,cont,SMF, ass_part, &
                            prod_boom, prod_rec,yi,yp,effi,effp,thetai,thetap,gend)
                                    
                f2  = valfres(1)
                
                 if(f1==-900000000000000000000000000d0 .AND. share == 0d0) write(*,*) 'vf is infinity (for couples) f1'
                
                
                d   =   alpha1*alpha2*d
                
              
            do while(d>tol)
                       
                d = d*alpha2

            if(f2<f1)then

                x2 = x1
                x1 = x1-d
                f2 = f1
                valfres  =  &
                valfun_sing(x1,share,gam,y,theta,eff,beta, &
                            c0,l,t,z,vfs,vc,pi_boom,pi_rec,marr,trans_part,nu,cont,SMF,ass_part, &
                            prod_boom,prod_rec,yi,yp,effi,effp,thetai,thetap,gend)
                                    
                f1 = valfres(1)
 
            else

                x1 = x2
                x2 = x2+d
                f1 = f2
                valfres  = &
                 valfun_sing(x2,share,gam,y,theta,eff,beta, &
                             c0,l,t,z,vfs,vc,pi_boom,pi_rec,marr,trans_part,nu,cont,SMF,ass_part, &
                             prod_boom,prod_rec,yi,yp,effi,effp,thetai,thetap,gend)
                                    
                f2 = valfres(1)
   
            end if
            
            end do    

            if(f2>f1) then
        
            x1=x2
            f1=f2
        
            end if
            
            ! store results in array
            output(1) = x1
            output(2) = f1
            output(3) = valfres(2)
            output(4) = valfres(3)
            output(5) = valfres(4)

    end function golden_sing
    
    ! 6) golden function (couple)

        function golden_coup(aa,bb,share,gami,gamp,yi,yp,thetai,thetap,effi,effp, &
                             betai,betap,c0,l,zi,zp,vfi, &
            vfp,vfc,vci,vcp,pii_boom,pii_rec,pip_boom,pip_rec,pic_boom,pic_rec,lambda,nu,cont,conti,contp,SMF, &
            yf,ym,efff,effm,thetaf,thetas,prodf_boom,prodf_rec, prodm_boom, prodm_rec) result(output_coup)
                                
       
            ! declare variables
            implicit none
    
            real*8  :: alpha1, alpha2, d, x1, x2, tol, output_coup(5), f1, f2
            
            real*8  :: aa, bb, r_intertmp(4)
            
            integer :: l,zi,zp
    
            real*8  :: share,c0,lambda
            
            real*8  :: prodf_boom(nzz),prodf_rec(nzz), prodm_boom(nzz),prodm_rec(nzz)
            
            real*8  :: yi, yp, thetai, thetap, effi, effp, gami, gamp, betai,betap, nu
    
            real*8  :: SMF, yf, ym, efff, effm, thetaf, thetas

            real*8, dimension(ncc,nzz) :: vfi, vfp
            
            real*8, dimension(ncc,nzz,nzz) :: vfc, vci, vcp 
            
            real*8, dimension(ncc) :: cont, conti, contp
            
            real*8, dimension(nzz,nzz) :: pii_boom, pii_rec, pip_boom, pip_rec
    
             real*8, dimension(nzz*nzz,nzz*nzz) :: pic_boom, pic_rec
                                        
            ! find optimal value with golden section search
    
            tol     =   0.000000001d0
            alpha1  =   (3d0-sqrt(5d0))/2d0
            alpha2  =   (sqrt(5d0)-1d0)/2d0
                d   =   bb-aa
                x1  =   aa+alpha1*d    
                x2  =   aa+alpha2*d    

                
                r_intertmp = valfun_coup(x1,share,gami,gamp,yi,yp,thetai,thetap,effi,effp, &
									     betai,betap,c0,l,zi,zp,vfi,vfp,vfc,vci,vcp, &
									     pii_boom,pii_rec,pip_boom,pip_rec,pic_boom,pic_rec, &
									     lambda,nu,cont,conti,contp,SMF,yf,ym,efff,effm, &
									     thetaf,thetas,prodf_boom, prodf_rec, prodm_boom, prodm_rec)
                        
                f1  =   r_intertmp(1)
                
                if(f1==-900000000000000000000000000d0 .AND. share == 0d0) write(*,*) 'vf is infinity (for couples) f1'

                r_intertmp = valfun_coup(x2,share,gami,gamp,yi,yp,thetai,thetap,effi,effp, &
									     betai,betap,c0,l,zi,zp,vfi,vfp,vfc,vci,vcp, &
									     pii_boom,pii_rec,pip_boom,pip_rec,pic_boom,pic_rec, &
									     lambda,nu,cont,conti,contp,SMF,yf,ym,efff,effm, &
									     thetaf,thetas,prodf_boom, prodf_rec, prodm_boom, prodm_rec)
                        
                f2  =   r_intertmp(1)
                
                if(f2==-900000000000000000000000000d0 .AND. share == 0d0) write(*,*) 'vf is infinity (for couples) f2'
                
                d   =   alpha1*alpha2*d

            do while(d>tol)
        
                d = d*alpha2
        
            if(f2<f1)then
        
                x2 = x1
                x1 = x1-d
                f2 = f1
                r_intertmp = valfun_coup(x1,share,gami,gamp,yi,yp,thetai,thetap,effi,effp, &
									     betai,betap,c0,l,zi,zp,vfi,vfp,vfc,vci,vcp, &
									     pii_boom,pii_rec,pip_boom,pip_rec,pic_boom,pic_rec, &
									     lambda,nu,cont,conti,contp,SMF,yf,ym,efff,effm, &
									     thetaf,thetas,prodf_boom, prodf_rec, prodm_boom, prodm_rec)
                        
                f1  =   r_intertmp(1)
            else
        
                x1 = x2
                x2 = x2+d
                f1 = f2
                r_intertmp = valfun_coup(x2,share,gami,gamp,yi,yp,thetai,thetap,effi,effp, &
									     betai,betap,c0,l,zi,zp,vfi,vfp,vfc,vci,vcp, &
									     pii_boom,pii_rec,pip_boom,pip_rec,pic_boom,pic_rec, &
									     lambda,nu,cont,conti,contp,SMF,yf,ym,efff,effm, &
									     thetaf,thetas,prodf_boom, prodf_rec, prodm_boom, prodm_rec)
                        
                f2  =   r_intertmp(1)
                
            end if
            
            end do    

            if(f2>f1) then
        
            x1=x2
            f1=f2
        
            end if
    
            
            ! store results in array
            r_intertmp=valfun_coup(x1,share,gami,gamp,yi,yp,thetai,thetap,effi,effp, &
									     betai,betap,c0,l,zi,zp,vfi,vfp,vfc,vci,vcp, &
									     pii_boom,pii_rec,pip_boom,pip_rec,pic_boom,pic_rec, &
									     lambda,nu,cont,conti,contp,SMF,yf,ym,efff,effm, &
									     thetaf,thetas,prodf_boom, prodf_rec, prodm_boom, prodm_rec)
            
            
            output_coup(1) = x1
            output_coup(2) = r_intertmp(1)
            output_coup(3) = r_intertmp(2) 
            output_coup(4) = r_intertmp(3) 
            output_coup(5) = r_intertmp(4) 


        end function golden_coup
        


end module
